Reassessing the formation of the inner Oort cloud in an embedded star cluster 
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r"| We re-examine the formation of the inner Oort comet cloud while the Sun was in its birth cluster with the aid of numerical 
Q-i simulations. This work is a continuation of an earlier study (Brasser et al., 2006) with several substantial modifications. First, the 
q system consisting of stars, planets and comets is treated self-consistently in our N-body simulations, rather than approximating the 
stellar encounters with the outer Solar System as hyperbolic fly-bys. Second, we have included the expulsion of the cluster gas, a 
c/5 feature that was absent previously. Third, we have used several models for the initial conditions and density profile of the cluster - 
i ^ i either a Hernquist or Plummer potential - and chose other parameters based on the latest observations of embedded clusters from 
the literature. These other parameters result in the stars being on radial orbits and the cluster collapses. Similar to previous studies, 
in our simulations the inner Oort cloud is formed from comets being scattered by Jupiter and Saturn and having their pericentres 
decoupled from the planets by perturbations from the cluster gas and other stars. We find that all inner Oort clouds formed in 
these clusters have an inner edge ranging from 100 AU to a few hundred AU, and an outer edge at over 100000 AU, with little 
variation in these values for all clusters. All inner Oort clouds formed are consistent with the existence of (90377) Sedna, an inner 
Oort cloud dwarf planetoid, at the inner edge of the cloud: Sedna tends to be at the innermost 2% for Plummer models, while it 
( ^ Js 5% for Hernquist models. We emphasise that the existence of Sedna is a generic outcome. We define a 'concentration radius' 
_ for the inner Oort cloud and find that its value increases with increasing number of stars in the cluster, ranging from 600 AU to 
t-H [1500 AU for Hernquist clusters and from 1500 AU to 4000 AU for Plummer clusters. The increasing trend implies that small star 
clusters form more compact inner Oort clouds than large clusters. We are unable to constrain the number of stars that resided in 
the cluster since most clusters yield inner Oort clouds that could be compatible with the current structure of the outer Solar System. 
J^~j The typical formation efficiency of the inner Oort cloud is 1.5%, significantly lower than previous estimates. We attribute this to 
rS jthe more violent dynamics that the Sun experiences as it rushes through the centre of the cluster during the latter's initial phase of 
violent relaxation. 
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1. Introduction and background 

There have been several studies of Oort cloud formation in a 
star cluster environment. In 1999 Soenke Eggers published his 
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Ph.D. thesis in which he analysed the effects of a star cluster on 
the formation and evolution of the Oort cloud (OC). He used a 
Monte Carlo method and two star clusters, in which the stellar 
encounters occurred at constant time intervals and their effects 
on the comets were computed analytically. The first cluster had 
an effective number density of 625 stars pc~ 3 and the other had 
an effective number density two orders of magnitude lower. For 
both clusters Eggers assumed a velocity dispersion of 1 km s _I . 
His model did not include a tidal field caused by the cluster 
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potential. Eggers denned a comet to be in the Oort cloud if 
it attained q > 33 AU and simultaneously a > 110 AU. With 
these definitions, he obtained efficiencies of 1.7% and 4.8% for 
the low and high density clusters respectively, with comets on 
orbits with a typical semi-major axis between 3 000 AU and 
6 000 AU. The clouds were also found to be mostly isotropic. 

A parallel study of the formation of the inner Oort cloud 
in a denser, clusteresque environment has been performed by 
Fernandez & Brumni (2000). They placed comets on eccentric 
orbits with semi-major axes ~ 200 AU and included an approx- 
imate model of the tidal field of the gas and passing stars from 
the cluster. The cluster had a maximum stellar number den- 
sity of 100 pc~ 3 and the maximum mass density of the core 
of the molecular cloud gas was 5000 M pc~ 3 . Their simula- 
tions formed a dense inner Oort cloud where the comets had 
semi-major axes of a few hundred to a few thousand astronom- 
ical units. The outer edge of this cloud was dependent on the 
density of gas and stars in the cluster. The most interesting 
part of their study was their ability to successfully deposit a 
fair amount of material that was scattered by Jupiter and Sat- 
urn, which were the main contributors to forming the inner Oort 
cloud. In the current environment, on the other hand, the con- 
tribution to the Oort cloud from Jupiter and Saturn is lower than 
that from Uranus and Neptune (e.g. Dones et al., 2004). How- 
ever, Fernandez & Brumni (2000) pointed out that if the Sun 
remained in this dense environment for long, the passing stars 
could strip a significant fraction of the comets away from the 
Sun and the trapping efficiency might end up being low. This 
low trapping efficiency was partially a result of their cluster life- 
times being too long. 

The interest in the formation of the inner Oort cloud in a 
cluster environment gained renewed interest with the discovery 
of (90377) Sedna, a dwarf planet with semi-major axis 500 AU 
and perihelion of 76 AU, so that its orbit is detached from that of 
Neptune (Brown et al., 2004). Gladman et al. (2002) has shown 
that the object 2000 CR105, having an orbit with q = 44 AU and 
a ~ 200 AU, could not be reproduced via chaotic diffusion, 
which ceases beyond q — 38 AU. Thus another mechanism 
had to be responsible for placing both 2000 CR105 and Sedna 
on their current orbits. Morbidelli & Levison (2004), together 
with Kenyon & Bromley (2004), successfully demonstrated 
that the most viable way to reproduce the orbits of Sedna and 
2000 CR105 is via a slow, close passage of a relatively heavy 
star. Morbidelli & Levison (2004) argued that the encounter 
had to happen early in order to still have a reasonably popu- 
lated Oort cloud afterwards. However, the low velocity of the 
encounter is difficult to obtain in the current Galactic environ- 
ment, and they suggested that this passage occurred while the 
Sun was in its birth cluster. 

The above results led to Brasser et al. (2006) - henceforth 
BDL6 - to investigate the formation of the inner Oort cloud 
in an embedded cluster environment, in which the gas from 
the molecular cloud is still present (Lada & Lada, 2003). In- 
spired by Fernandez & Brumni (2000) they attempted to con- 



strain the environment that was needed to save comets under 
the dynamical control of Jupiter. They employed a Plummer 
model (Plummer, 1911) to construct a series of clusters with 
varying central density but with a more or less fixed number 
of stars, because most stars form in clusters of a few hundred 
stars (Lada & Lada, 2003). BDL6 used a simple leapfrog inte- 
grator for the cluster in the Plummer potential and recorded the 
positions and velocities of stars, including time, as they came 
within a user-specified distance of the Sun. The encounter data 
were then used in SWIFT RMVS3 (Levison & Duncan, 1994). 
The latter was modified to include the effects of the stars, as in 
Dones et al. (2004), and the gravitational force from the cluster 
gas on the comets and the planets. The Sun was assumed to be 
on a fixed orbit in the Plummer potential from which the tidal 
torque of the gas on the comets was computed. The typical effi- 
ciency for the formation of the inner Oort cloud was 10%. They 
showed that Sedna's orbit could be reproduced when the cen- 
tral density of the cluster exceeded 10000 M Q pc~ 3 in gas and 
stars. For these clusters Sedna was found to be at the inner edge 
of the inner Oort cloud. In order to reproduce the orbit of 2000 
CR105, an even higher (central) density was needed. However, 
the orbital distribution of the inner Oort cloud in these very high 
density clusters is found to be inconsistent with the current ob- 
servations of the outer solar system (Schwamb et al., 2010), so 
that the clusters where Sedna is at the inner edge are preferred. 

In a similar study, Kaib & Quinn (2008) studied inner Oort 
cloud formation in an open cluster with stellar number densities 
ranging from 10 pc~ 3 to 100 pc~ 3 . Kaib & Quinn (2008) mod- 
elled the effect of the stars using the same approach as Dones et 
al. (2004) and BDL6. The maximum cluster life time was set 
to 100 Myr and the density of the stars in the cluster decayed 
linear with time, to account for mass loss by mutual scattering 
of the stars. When the cluster had completely disappeared Kaib 
& Quinn (2008) continued their Oort cloud simulations until the 
age of the Solar System, something BDL6 did not do. Quanti- 
tatively their results were similar to BDL6 and they were able 
to reproduce Sedna when the stellar number density exceeded 
30 pc- 3 . 

In a recent publication attempting to solve some of the out- 
standing problems associated with the Oort cloud as a whole, 
Levison et al. (2010) investigated the capture by the Sun of 
comets from other stars. They simulated embedded clusters 
ranging from 30 to 300 stars with a star formation efficiency 
of 10% to 30%. They placed a disc of comets around each star 
with random orientation and orbits with q — 30 AU and semi- 
major axes ranging from 1 000 to 5 000 AU. The whole sys- 
tem of stars and comets was simulated until the median spacing 
between the stars became 500000 AU. From these numerical 
simulations Levison et al. (2010) concluded that the capture ef- 
ficiency is high enough to obtain the current population of the 
Oort cloud provided that most of the stars in the cluster con- 
tained a similar number of comets to the Sun. At least 90% of 
the comets in the Oort cloud could be extrasolar in origin. 

Unfortunately, apart from the Levison et al. (2010) study, 
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which relied on extrasolar comets rather than indigenous comets 
to populate the Oort cloud, all of the above works suffer from 
the limitation that the stars are treated as hyperbolic encoun- 
ters and the Sun is on a fixed orbit in the cluster. Both of these 
assumptions are wrong: the stars' motion with respect to the 
Sun reverse direction when their distance to the Sun is of the 
same order as the size of the current Oort cloud. If the stellar 
number density in a cloud is n* pc~ 3 , then their average nearest- 
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neighbour distance is r — 0.62n„ ' pc, and taking a typical 
density of n, = 30 pc" 3 yields r = 0.2 pc or 42 000 AU. Sec- 
ondly, mutual scattering among the stars changes their orbits 
and some end up on highly elliptical orbits on their way to be- 
ing ejected. Hence the assumption of a static orbit is no longer 
valid. A third issue is gas removal. BDL6 stopped their sim- 
ulations after 3 Myr by assuming that at this time the Sun left 
the cluster. Kaib & Quinn (2008) decreased the density of their 
fictitious open cluster linearly and it was gone after 100 Myr. 
Levison et al. (2010) made the gas go away exponentially with 
an e-folding time of 10000 yr. Fourth, the BDL6 study relied 
on very high central gas densities in order to torque comets un- 
der the dynamical control of Jupiter into the inner Oort cloud 
before they were ejected. While BDL6 argued that the densi- 
ties that were chosen are in agreement with the peak densities 
observed in some embedded clusters (Gutermuth et al., 2005), 
their initial conditions can be improved by using more recent 
observational data and better models for the embedded clusters. 
Any reasonable model of the formation of the inner Oort cloud 
in a star cluster environment has to take the above issues into 
account. The aim of this paper, therefore, is to re-investigate the 
formation of the inner Oort cloud in a cluster environment us- 
ing (i) a better model for the embedded star cluster which best 
matches the current observations, in particular the initial condi- 
tions for the stars, gas and the gas dispersal, and (ii) a computer 
code that can handle stars, planets and comets at the same time 
so that there is no need to rely on the assumption of hyperbolic 
stellar encounters. 

Levison et al. (2010) used a computer code, based on SyMBA 
(Duncan et al., 1998), that was able to integrate both the comets 
and stars symplectically and self-consistently without relying 
on assumptions of hyperbolic flybys. In this study we shall 
use their code. This paper is divided as follows. In Section 2 
we summarise some of the basic properties of embedded star 
clusters that we need for our simulations as inferred from ob- 
servations. Section 3 deals with the initial conditions and meth- 
ods of our numerical simulations. In Section 4 we present the 
properties of the inner Oort cloud resulting from the numerical 
simulations. In Section 5 we compare these results with recent 
observations of the outer Solar System. In Section 6 we present 
our discussion, followed by the conclusions in Section 7. 

2. Cluster properties and models 

In this section we discuss the models and parameters that 
we employed for the simulation of the embedded star clusters. 
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Table 1: Clusters common to Gutermuth et al. (2009) and Lada & Lada (2003) 
for which the radius is known. The columns are: name, total radius from Guter- 
muth et al. (2009), core radius from Gutermuth et al. (2009), radius from Lada 
& Lada (2003), total number of stars from Gutermuth et al. (2009), number 
of stars in the core from Gutermuth et al. (2009) and number of stars listed in 
Lada & Lada (2003). 

2.1. Cluster size 

We use the data of embedded star clusters within 2 kpc of 
the Sun from Lada & Lada (2003) and Gutermuth et al. (2009). 
Adams et al. (2006) use the Lada & Lada (2003) data and find 
that most embedded clusters have between 100 to 1 000 mem- 
bers. The cumulative distribution of the number of stars, N, has 
a best fit / = 0.637 log/Y - 1.045 i.e the distribution increases 
logarithmically with N (Adams, 2010). 

For most embedded clusters, the observed surface density 
is a constant (Allen et al., 2007) i.e. N/R 2 is constant, where 
R is the radius of the cluster. This relation is in agreement 
with the observed density structure of giant molecular clouds 
(Blitz et al., 2007) and theoretical modelling of cloud collapse 
(Larson, 1985). The observed constant surface density implies 
that the size of the cluster scales with the number of stars as 
R = R N 1/2 , where R Q is a scaling parameter. From the cata- 
logue of Lada & Lada (2003), Adams et al. (2006) find a best 
fit R = R (N/I00) l/2 pc where R e (0.577, 1) pc. For their 
subsequent simulations Adams et al. (2006) use Rq = 0.577 pc 
to maximise dynamical interactions among the stars. However, 
the radii of embedded clusters presented in Lada & Lada (2003) 
might only be valid for the cores of the clusters that are listed. 
Recently Gutermuth et al. (2009) performed a Spitzer study of 
a large sample of embedded clusters, some of which are also 
listed in Lada & Lada (2003). Gutermuth et al. (2009) char- 
acterise each cluster by a core and an extended 'halo'. Us- 
ing nearest-neighbour distance counts and assigning a radius 
to the cluster as being half the distance between the two far- 
thest stars, a best fit through their data for the halos yields 
R H = Ro(N/lOOf where R = 1.92 ± 0.52 pc and the expo- 
nent fi = 0.41 ± 0.09. For the cores R c = R (N /I00f with 
R = 0.95 ± 0.36 pc and /3 = 0.47 ± 0.11. The fitting param- 
eters and their error bars are displayed in Fig.Q] The common 
clusters of Gutermuth et al. (2009) and Lada & Lada (2003) for 
which radii are available are listed in TableQ] The columns are: 
name, total radius from Gutermuth et al. (2009), core radius 
from Gutermuth et al. (2009), radius from Lada & Lada (2003), 
total number of stars from Gutermuth et al. (2009), number of 
stars in the core from Gutermuth et al. (2009) and number of 
stars listed in Lada & Lada (2003). Clusters which are common 
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Figure 1: Range in sizes, R, vs exponents, /?, for the cluster core sizes from 
Lada & Lada (2003) and Gutemiuth et al. (2009). 



to both catalogues but have no radius listed in Lada & Lada 
(2003) are AFGL 490, IC 5136, Cep C, LH Ha 101, Serpens, 
Cep A, L988-e and R CrA. As can be seen from Table Q] most 
of the sizes listed in Lada & Lada (2003) are close to the size 
of the cores listed in Gutermuth et al. (2009), or are interme- 
diate between the core and the halo. In any case, the cluster 
sizes in Lada & Lada (2003) are systematically smaller than in 
Gutermuth et al. (2009), and the fit through the core data of 
Gutermuth et al. (2009) is compatible with the fit through the 
data of Lada & Lada (2003), but the halos are not. 

The best fits seem to indicate that R ~ N p with >S e [1/3, 1/2]. 
However, the A^ 2 relation is based on the assumption that the 
column or surface density of stars is more or less constant and is 
an artefact of the way the stars are counted. Most star identifica- 
tion algorithms rely on the density-weighted nearest-neighbour 
method of Casertano & Hut (1985), or the minimum spanning 
tree method (e.g. Graham & Hell, 1985). Both are often em- 
ployed to identify star clusters (Bastian et al., 2007, 2009; Cartwri| 
& Whitworth, 2004; Gutermuth et al., 2009; Schmeja & Klessen, 
2006). However, in all cases the cluster radius is defined as the 
radius of a circle with the same area as the projected cluster 
(Schmeja, 2011). All methods truncate the size of the cluster 
when the projected distance between two neighbouring stars is 
larger than some threshold value, which is equivalent to assum- 
ing that the surface density inside said circle is more or less 
constant. Thus, the size of the cluster then obviously scales 
as A^ 2 . However, the N 1 ^ 2 relation appears in disagreement 
with another observation, and that is that the average stellar 
number density in these clusters is more or less constant (Car- 
penter, 2000; Lada & Lada, 2003; Proszkow & Adams, 2009), 
with a median number density of Hm = 65 pc 3 . From the 
clusters listed in Lada & Lada (2003) and Gutermuth et al. 
(2009), we compute the median value for the halos of Guter- 
muth et al. (2009) to be um = 3.1 pc 3 while for the cores it is 
«m = 46.2 pc 3 , comparable to the value listed earlier. How- 
ever, the average stellar number density from one cluster to the 
next can vary by approximately an order of magnitude. This 



true for both cores and halos. Thus the values quoted above 
should be interpreted as indicative only. 

It is easy to verify whether the observation of more or less 
constant stellar number density is consistent with the relation 
between the cluster's size and the number of stars. The total 
number of stars in the cluster, apart from a constant, is N - 
(n)R 3 , where (n) is the average stellar number density. Substi- 
tuting R — N@ we have (n) = N l ~ 3/3 , which is a constant only 
if B = 1/3. Thus it seems sensible to adopt a cluster size that 
scales as N l/3 . In order to determine whether or not this scaling 
makes sense, we turn to observations of open and globular star 
clusters for guidance. 

King (1962) proposed that the size of a star cluster, whether 
open or globular, is given by its tidal radius, r,. At this distance 
from the centre, the tidal effects of the Milky Way Galaxy start 
to dominate over the self gravity of the cluster. King (1962) 
states that 



n \4A(A-B)/ ' 



(1) 



where Q is the gravitational constant, M is the total mass of 
the cluster and A and B are the Oort constants (e.g. Binney 
& Tremaine, 1987). The tidal radius scales as N 1 ^ 3 because 
M oc N. When adopting a flat Galactic rotation curve with 
angular velocity Qc = 30 km s _1 kpc~' (MacMillan & Bin- 
ney, 2010), we have A — \B\ — 15 km s kpc~' and so r, ~ 
4.6 (N/ 100) 1/3 pc, which is much larger than the halo sizes 
of Lada & Lada (2003) and Gutermuth et al. (2009). This 
discrepancy is most likely caused by the fact that for embed- 
ded clusters the background density is not that of the Galac- 
tic disc, but rather of the surrounding molecular cloud. Typ- 
ical densities of molecular clouds are some ~ 1 M pc~ 3 , so 
that the tidal radius in equation ([TJ above should be divided by 
~ 2 i.e. r, ~ 2.3 (A^/100) 1 ^ 3 pc. This result is similar to and 
compatible with the halo distances obtained by Gutermuth et 
al. (2009). The shape of the zero-velocity curves of the clus- 
ter reduce the tidal radius even further (Innanen et al., 1983) 
to r, ~ 1.7 (A7100) 1 '' 3 pc. Given the uncertainties in the ob- 
served properties of the clusters and in the size versus num- 
ber of stars, we shall anchor the value of r, for N = 100 to 
the best-fit halo value of Gutermuth et al. (2009) and thus use 
r, = 1 .92 (A7 100) 1/3 pc for the size of the cluster. 



2.2. Core radius and internal structure 

Open and globular star clusters usually show a dense core 
with more or less constant surface brightness and then an ex- 
tended halo where the surface density falls off (King, 1962), 
typically as r~ r , where y e (2,3) (King, 1962; Elson et al., 
1987). Observers generally define the cluster core radius, r c , as 
the distance from the cluster centre at which the surface bright- 
ness drops by a factor of two from the central value. Theo- 
rists, however, often use a definition based on the central den- 
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sity, po, and central velocity dispersion, cr . The core radius is 
then given by (King, 1966) 

2 V^Gpo ' *" ^ 

which we shall use here. For most cluster models the core ra- 
dius corresponds roughly to where the stellar volume density 
has decreased by a factor of 3. A second method computes a lo- 
cal density using a star's nearest neighbours (Casertano & Hut 
1985), in which the core radius becomes a density-weighted 
quantity obtained from the root-mean-square stellar distances. 
We refer to Portegies Zwart et al. (2010) for a more in-depth 
discussion on how the core radius is defined and measured. 
King (1962) defines the 'concentration ratio' of the cluster as 
c K = \og(r,/r c ), where the log is with base ten; we shall use 
the non-logarithmic form c = r t jr c here. For young open clus- 
ters (Piskunov et al., 2008) the value of c is typically 3 to 6. 
The embedded cluster table of Gutermuth et al. (2009) yields 
similar values. For the young open, non-relaxed cluster NGC 
6611, whose estimated age is 1.3 Myr, c ~ 10 (Bonatto et al., 
2006). All of these observations suggest that c ranges from ap- 
proximately 3 to 10, and thus these clusters have fairly shallow 
profiles and central potential wells. We shall use c — 3 and 
c = 6 in our simulations. 



Adams et al. (2006) set a H = R, with R the size of the clusters 
obtained from Lada & Lada (2003). The total mass can be con- 
verted to a 'central density' through po = MI{2na i H ). The total 
mass inside a H is jM and here p(r) ~ r~ x . In order to model 
clusters where p(r) ~ r° close to the centre, one could use a 
Dehnen profile with y — 0. However, we decided to settle for 
the density profile of Plummer (1911), which is widely used in 
star cluster simulations because of its simplicity (e.g. Aarseth 
et al. 1974; Baumgardt & Kroupa, 2007; Kroupa et al., 2001). 
Its density profile is given by 



Pp(r) 



3M 



-5/2 



(6) 



where the central density is po = 3M/(4-na 3 p ) and the density at 
the centre scales as p(r) ~ r°. Here a P is the Plummer radius. 
The potential is given by 



. . GMI r 2 \~m 

<r) = 1 + -=- 

flp V ail 



(7) 



Here we shall use both the Hernquist and Plummer distributions 
only for their simplicity and ability to adequately reproduce the 
observed density structure in the centre of the cluster. 



Inside the cluster the volume density scales as p(r) oc r~ y 
where y e (0,2) (Schmeja & Klessen, 2006; Schmeja et al., 
2008; Andre et al., 2007), with values in the range to 1 being 
the most common. It is well known that the density of most 
globular clusters are best fitted with King profiles (King, 1966), 
which have a well-defined core and halo. Inside the core the 
density is more or less constant while outside the core the den- 
sity falls off quickly (Portegies Zwart et al., 2010). However, it 
is unclear if the King profiles are suitable for young/embedded 
clusters (Portegies Zwart et al., 2010). Since the potential for 
the King models cannot be written in closed form as a function 
of distance from the centre, which we need in our computer 
code, we prefer not to use these models. There exist various 
alternatives in the literature to compute the density and poten- 
tial. For spherical Galaxies the models by Dehnen (1993) and 
Tremaine et al. (1994) are often used. The volume density of 
these profiles are 

(3 - y)M 



P(r) = 



a D 



(3) 



4n r?(r + a ) 4 ~ y ' 
where M is the total mass in gas and stars, ap is a parameter 
radius and y measures the density concentration at the centre. 
The density profiles of Jaffe (1983) and Hernquist (1990) are 
the cases with y — 2 and y = 1 respectively. In their clus- 
ter simulations Adams et al. (2006) use the density profile of 
Hernquist (1990), for which the density is given by 

M a H 1 



(4) 



2n r (r + afj) 3 ' 
where a H is the Hernquist radius. The corresponding potential 



is 
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An additional quantity to address is the magnitude of the 
velocity dispersion within the clusters. Observations indicate 
that in the youngest embedded clusters the velocities of star- 
less clumps and young stellar objects are a fraction of the virial 
value. Thus, the orbits of stars are mostly radial. Hydrody- 
namical simulations of cluster formation from dynamically hot 
gas results in the formation of young stellar objects that have 
speeds comparable to the sound speed, which are much lower 
than the virial value (Bate et al., 2003). An example of a cluster 
with sub-virial speeds is L1688, part of the p Oph complex, 
where the velocities of the stars are approximately 30% of the 
virial value (Andre et al., 2007). Similar results are found in p 
Oph A (Di Francesco et al., 2004) at 50% of the virial value, 
while NGC 2264 (Peretto et al., 2006) and NGC 1333 (Walsh 
et al., 2007) have even lower values. Indeed, some of the pre- 
stellar cumps and some young stellar objects in these clusters 
appear to exhibit a collapse because of the low speeds. The 
deviation from virial equilibrium of a cluster is quantified by 
the parameter Q, which is the ratio between the total kinetic en- 
ergy and total potential energy; virial systems have Q-\. For 
our purpose we shall consider a velocity dispersion with a value 
between 0.3 and 0.5 of the virial value i.e. Q e [0.05,0.125]. 

The last issue we discuss is mass segregation. Massive stars 
are believed to sink to the cluster center over a relatively long 
time scale, given approximately by ?r/M*, where f« is the dy- 
namical relaxation time, and M, is the mass of the star in solar 
masses (Portegies Zwart, 2009). Of course, the massive stars 
can also be formed at the cluster centres, and some observa- 
tional evidence (Testi et al. 2000; Peretto et al., 2006) and theo- 
retical considerations (Bonnell & Davies 1998; McKee & Tan, 
2003) support this point of view. However, Allison et al. (2009) 



report very rapid mass segregation if the initial structure of the 
cluster is sufficiently fractal and the velocities are highly sub- 
virial, typical for these young clusters. In their study of these 
highly fractal clusters with subvirial velocities, mass segrega- 
tion was achieved for stars heavier than 5 M Q and completed in 
approximately 1 Myr. This rapid mass segregation appears con- 
sistent with observations of the Orion Nebula Cluster (Hillen- 
brand & Hartmann, 1998; Moeckel & Bonnell, 2009). In this 
study we consider clusters which are already mass segregated 
for stars heavier than 5 M Q . It should be noted that the rapid 
mass segregation coincides with a violent relaxation phase af- 
ter the gravitational collapse, both of which are a result of its 
subvirial velocities, as the cluster tries to reach equipartition. 
As it does so it shrinks in size by approximately a factor two or 
more (Allison et al., 2009). In this study we do not concern our- 
selves with the mechanism behind mass segregation, whether it 
is through dynamics or by formation, but assume it has already 
happened for stars heavier than 5 M e . 

Now that we have discussed most of the properties of em- 
bedded clusters based on the latest observations, and constrained 
some of the key parameters that will be used, we next describe 
our numerical methods. 

3. Initial conditions and numerical methods 

In this section we describe the initial conditions and meth- 
ods employed for our numerical simulations. 

3.1. Initial conditions and gas removal 

We generate the stars in each cluster as follows. First, the 
desired number of stars was chosen. The mass of each star was 
then calculated randomly according to the Initial Mass Function 
formulation of Kroupa et al. (1993), with the functional form 



M(£) = 0.08 + 



0.19£ L55 + 0.05£° 



).6 



(W) 1 



,0.58 



■ Mp 



(8) 



Here £ is a number chosen randomly on the interval [0, 1) and 
M(£) is the mass of the star in solar masses. The average stellar 
mass is then <m»> = / M(£)d£, ~ 0.43 M Q . The total mass of 
the stars is approximately M, = {m*)N, and is valid for large 
N. No primordial binaries were included. In order to model 
the Solar System we generated one star with a mass of exactly 
1 Mq and considered it to be the Sun. The tidal radius of the 
cluster was computed as r, = 1.92 (N/ 100) 1/3 pc, and the core 
radius is either or |r f . For the Plummer profile the Plummer 

radius, a P , is then computed from r c by using cr = 1 and 

solving for ap, resulting in ap — V2r c . For the Hernquist pro- 
file, the central velocity dispersion is 0, so we use p(r c ) = ^po 
and solve for a H to find a H « 1.46r e . 

The stars in the cluster are subjected to three forces. The 
first is their mutual gravitational interaction. The second is 
caused by the Galactic tide and bulge, and is modelled accord- 
ing to the formulation of Levison et al. (2001) but with the 



Oort constants set at A = \B\ = 15 km s _1 kpc~' respectively 
(MacMillan & Binney, 2010). The third force is caused by the 
gas that is present in the cluster, whose density profile is also 
modelled either by the Hernquist or Plummer distribution, with 
the same values of a H or a P used for the stars and gas. The 
mass in gas is related to the total mass in stars, M„ and the star 
formation efficiency (sfe, s) by M g = (e _1 - \)M„. The sfe was 
set either to 0.1 or to 0.25, which mostly brackets the observed 
range (Lada & Lada, 2003). 

The magnitude of the position and velocity vectors of the 
stars are generated from the isotropic energy distribution func- 
tions of either the Hernquist or Plummer profiles using a von 
Neumann rejection technique (Press et al., 1992). For the Plum- 
mer model we followed Aarseth et al. (1974). They solve the 
equation M(r) = £M for r for each star, with f = nM*/N 
and n < N. This results in successive values of M(r) being 
evenly spaced, which appears to be fine for the heavy stars that 
have settled in the centre, but is artificial for the other stars. 
Thus we decided to adopt this method for the heavy, segregated 
stars but replaced £ by a random number on the interval be- 
tween [0,£ max ) for the other stars. Here £ max is determined by 
M(r) = M(r t ). For the Hernquist model we followed the same 
procedure. The singularity of the distribution function of the 
Hernquist model at energy E = -GM/an is avoided by noting 
that 4nr 2 v 2 f(\E\) < 1.1, where f(E) is the Hernquist distribu- 
tion function and v is the velocity of a star. We use the latter 
formulation with the von Neumann technique (D. C. Heggie, 
personal communication). The position and velocity vectors 
were calculated with random orientation. In order to avoid the 
Sun leaving the cluster immediately, we 'rigged' the system by 
requiring that the Sun moves inwards if it is farther then jr, 
from the centre when the simulation is started. This was ac- 
complished by requiring that for the Sun r • v < if r > \r t . 

After generating the velocities and positions from the dis- 
tribution function the kinetic energy of the stars was reduced. 
A single value for the virial parameter Q for each simulation 
was randomly computed on the interval Q e [0.05, 0.125], and 
the kinetic energy of each star was then reduced accordingly. 
The maximum speed of the stars is v = (20'' 2 v esc where v esc = 
|20(r)| 1/2 is the local escape speed and O(r) is the gravitational 
potential of the Hernquist or Plummer sphere. The stars' low 
kinetic energy implies that the stars are initially on nearly ra- 
dial orbits and exhibit almost a free fall (which justifies having 
the Sun move inwards at first). The resulting set of stellar posi- 
tions and velocities was then used as the starting conditions in 
our full simulations. 

In the above procedure, all stars were assumed to have formed 
at the same time. There is an ongoing debate about whether or 
not this is true, but based on observations of NGC 2264 Peretto 
et al. (2006) favour the turbulent formation mechanism of Mc- 
Kee & Tan (2003), in which most stars form within ~ 10 5 yr, 
with the heaviest stars in the centre before the initial collapse of 
the cluster. Palla & Stahler (1999) report a similar conclusion 
based on the Orion cluster while Murray (201 1) states that most 
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stars form within the free-fall time of the cluster. These results 
suggest that most stars in the cluster form within a short time of 
each other, justifying our procedure. 

The next thing to model is the decay of the gas. The early 
work by Lada et al. (1984) removed the gas either instanta- 
neously or on a time scale ranging from three to four crossing 
times. The crossing time is defined as the time it takes for 
a star to cross the whole cluster i.e. t c = 2R/cr, where <x is 
the velocity dispersion. For our clusters t c ~ 4 Myr. Adams 
et al. (2006) kept the gas density constant and then removed 
it instantaneously after 5 Myr, independent of the properties 
of the cluster. Baumgardt & Kroupa (2007) removed the gas 
on a time scale from zero to ten crossing times. Proszkow & 
Adams (2009) removed the gas instantaneously at times rang- 
ing from 1 Myr to 7 Myr, to account for the spread in observed 
embedded cluster lifetimes. Levison et al. (2010) kept the gas 
density constant for 3 Myr and removed it exponentially with 
an e-folding time of 10000 yr. Thus, previous studies show 
great variability in their choice of the time of gas removal and 
its decay rate. An initially bound cluster responds to external 
perturbations and changes in the potential in approximately a 
crossing time and thus this time scale serves as a bench mark. 
Kroupa (2000) reports that typically the gas in an embedded 
cluster is removed on a crossing time. However, for our pur- 
poses we do not want to end up with a bound system of which 
the Sun is a potential member, and thus we would like to take 
the gas away quickly to ensure that the Sun escapes from the 
system. Even though Proszkow & Adams (2009) took the gas 
away instantaneously, they frequently found that bound systems 
remained. The strongest dependency appears to be on the star 
formation efficiency, s, but even when setting s = 0.1, they re- 
ported that some 15% of systems remained bound after the gas 
went away. Ultimately, the gas in the cluster is removed be- 
cause stellar winds from heavy stars create Ha bubbles and su- 
pernovae. Stellar winds cause rapid outflow of Ha bubbles, typ- 
ically with velocities of v w ~ 25 km s _1 (Whitmore et al., 1999). 
This leads to a removal time scale of td — R/v w ~ 80 000 yr, for 
clusters with N e [10, 1000]. Given this rapid time scale, we 
have decided to proceed as follows: the gas is kept at its ini- 
tial density for a time of either 2 Myr or 4 Myr. Given that the 
typical lifetime of an embedded cluster is some 5 Myr (Lada & 
Lada, 2003) with a maximum of 10 Myr, and that circumstel- 
lar discs have a median lifetime of 3 Myr (Currie et al., 2008), 
and that Jupiter and Saturn had to form in this time (Lissauer & 
Stevenson, 2007), a typical cluster lifetime after the formation 
of Jupiter and Saturn of 1 Myr to 5 Myr is reasonable. After the 
constant density phase, the density of the gas decays exponen- 
tially with an e-folding time of td = 30 000 yr. 

3.2. Numerical integration method 

Star clusters are not well-suited to the symplectic integra- 
tion methods typically used to follow the orbits of comets about 
stars (Spurzem et al., 2009). These methods split the Hamil- 
tonian into the sum of a part representing Keplerian motion 
around a fixed star and a part representing the non-Keplerian 



perturbations from other stars (Wisdom & Holman, 1991). How- 
ever, we find that the orbits of the stars and other bodies are 
efficiently computed using the standard Leap-Frog integrator, 
which is a second-order symplectic integrator where the Hamil- 
tonian is split into the sum of a part representing the kinetic en- 
ergy and one representing the gravitational potential energy. As 
described in Duncan et al. (1998), the key idea is to incorpo- 
rate a multiple time step symplectic method such that whenever 
two or more stars, or a comet and star(s), or comet and planet 
suffer close encounters, the time step for the relevant bodies is 
recursively divided to whatever level is required to resolve the 
encounter. In other words, the code splits the encounter into 
a series of 'shells' which Duncan et al. (1998) place apart as 
Ri/Ri+i = 3 2 ^ 3 and the time step is divided by 3 when crossing 
into the next shell. Since these requirements are already im- 
plemented in the SyMBA integration package (Duncan et al., 
1998), we modified it to use the Leap-Frog scheme. From now 
on this new code will be referred to as SyMBAC. For a descrip- 
tion of tests of the code without planets we refer to Levison et 
al. (2010). 

Including the Jovian planets orbiting the Sun in SyMBAC 
proved to be tricky. The Leap-Frog integration scheme is not 
very accurate for Kepler orbits unless one takes of the order 
of 1 000 steps per orbit (Kokubo et al., 1998). However, even 
with such a small time step it introduces a secular drift in the 
argument of pericentre, lo, of the planets (Kokubo et al., 1998). 
This drifting in u changes the secular properties of the planets 
and it is easy to grow the eccentricities of Jupiter and Saturn to 
crossing values. One solution to this problem is to invent a new 
symplectic method to handle both Kepler orbits and stellar en- 
counters, or to search for a set of parameters that keep the plan- 
ets stable. We chose to do the latter. It turned out that it is im- 
perative to keep Jupiter and Saturn on the same 'shell', and that 
the shells are scaled as = 3 and the time step between 

shells is divided by 6. The first 'shell' around the Sun is set to 
2916 AU and the time step is 100 yr, and is used for all of the 
calculations that follow. Since the typical spacing between the 
stars is much larger than 3 000 AU, encounters between stars are 
relatively rare and the code does not often need to resort to the 
recursive encounter routines apart from in the beginning when 
there are many comets close to the Sun. Hence the simulations 
we performed were finished in mere hours rather than days or 
weeks. As in BDL6 we do not include Uranus and Neptune, 
because their formation mechanism is still not well understood 
and the time scale of their formation is not well constrained 
(Goldreich et al., 2004). 

Comets were added only around the Sun as bodies with 
masses thirteen orders of magnitude smaller than the stars and 
planets. The stars and planets mutually interact, but the comets 
are only affected by the stars and the planets, and not by each 
other. The comets were started on orbits between 4.5 AU and 
12 AU with eccentricities and inclinations between and 0.01 
(radians). A total of 4 000 comets were used in each simulation. 
In order to hasten the evolution, we removed comets from the 
simulation when they were farther than 1.5r, from the cluster 
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centre and unbound. In order to keep the calculations traceable, 
and because we are only interested in the formation of the in- 
ner Oort cloud with indigenous comets, no comets were added 
around other stars, unlike in Levison et al. (2010). The simula- 
tions were stopped either after some maximum time (15 Myr) 
or when the Sun was farther than 1.5r, from the cluster centre 
after the gas began to evaporate. 

The last parameter to choose is the number of stars in the 
cluster. We have chosen seven individual values which bracket 
the approximate observed range of embedded clusters: 50, 100, 
250, 350, 500, 750 and 1 000. A total of 40 simulations were 
performed for each cluster, with 10 for each combination of the 
time the gas density remained constant, t g , and the concentra- 
tion c. Together with the two values of the sfe and the two po- 
tential/density pairs, this results in a total of 1 120 simulations, 
which were all performed on the CRIMSON Beowulf cluster at 
the Observatoire de la Cote d' Azur. 

3.3. Sample cluster evolution 

In this subsection we give an example of the evolution of 
one of the many star clusters that we simulated. The cluster 
is of the Hernquist type with N = 250, e = 0.1, t g = 4 Myr 
and the concentration parameter c = 6. Figures [2] and [3] show 
snapshots of the stars in the xy-plane at times indicated in the 
top-right corners of the panels. The Sun has been coloured in 
green to distinguish it from the other stars, and the size of the 

1/3 

dots correlates to the mass of the star as M, . 

The top-left panel of Fig. [2] depicts the initial positions of 
the stars, projected onto the xy-plane. These were generated 
according to the methods described above. Since the stars are 
all sub-virial, the cluster undergoes a state of collapse, through a 
process called 'violent relaxation' (Lynden-Bell, 1967). Some 
embedded clusters are observed to be in this state of collapse 
(Peretto et al., 2006; Andre et al., 2007). As one can see in the 
top-right panel of Fig. [2] the cluster has shrunk significantly and 
the density in the core is much higher than it was initially. Af- 
ter about 1 .4 Myr, the cluster has reached its smallest size and 
maximum core density. The cluster has a size of approximately 
1 pc. After this maximum compression of the core, the cluster 
expands again (bottom-right panel) and some of the stars are 
either unbound or on very elongated orbits. Continuing on to 
Fig. [3] we see that the cluster continues to expand and reaches 
some sort of a steady-state. At 4 Myr the gas goes away and in 
the bottom two panels of Fig. [3] we see the stars escaping from 
the system. This sequence of events, in which the cluster col- 
lapses and then expands again, occurs in all our simulations. 

Now that we have given a brief overview of the dynami- 
cal evolution of the star clusters, we turn to the results of our 
simulations of stars and comets. 

4. The formation of the inner Oort cloud 

In this section we present the results of our numerical simu- 
lations. Many properties of the inner Oort cloud that are formed 



in this way are similar to earlier results presented in Fernandez 
& Brunmi (2000), BDL6 and Kaib & Quinn (2008), such as 
the the cloud being nearly isotropic apart from the inner 10% 
or so. We shall not repeat all of these here but instead only 
focus on the key aspects and how these scale with the clus- 
ter properties: the distribution in semi-major axis, inclination 
and perihelion, the location of the inner edge and the formation 
efficiency. During this investigation it became clear that the pa- 
rameter space is very large, so we shall try to reduce it first. 
When examining the output from the simulations, it turned out 
that the initial virial ratio, Q, does not impact the results beyond 
the statistical noise one expects from one simulation to another. 
This result is surprising, because the value of Q determines how 
much the cluster shrinks during the initial collapse and violent 
relaxation. Thus, in what follows we consider the results to be 
averaged over the value of Q. This leaves us with the number 
of stars, N, the star-formation efficiency, e, the time the gas was 
removed, f^, the concentration, c, and the density profile (Plum- 
mer or Hernquist). For reference, a comet is defined to be in the 
cloud if it satisfies both a > 50 AU and q > 35 AU. 



4.1. The size and the concentration radius of the inner Oort 
cloud 

A simple way to depict the size and distance of the inner 
Oort cloud to the Sun is to plot the cumulative semi-major axis 
distribution of comets that are in the cloud. We have plotted 
this distribution for Hernquist clusters with 50, 250, 500 and 
1 000 stars (Fig. |4]l and Plummer (Fig. |5]l distributions respec- 
tively. Both figures correspond to clusters with an sfe of 10%. 
The red lines are for clusters with t g = 2 Myr and c = 3, green 
lines have t g = 4 Myr and c — 3, the blue lines correspond to 
the parameters t g = 2 Myr with c = 6, while the magenta lines 
depict cases with t g = 4 Myr and c = 6. We will use these 
t g - c-colour correlations throughout the paper, unless specified 
otherwise. The plots show a few interesting features that re- 
quire further discussion. 

The first is that, generally but not exclusively, the magenta 
lines precede the blue lines, which precede the green lines, 
which in turn precede the red lines. In other words, the inner 
Oort clouds formed above become more centrally condensed 
when both t g and c increase. This is easy to imagine. As the 
Sun resides in the cluster for longer times, the outermost comets 
are stripped away by encounters with the other stars. Only the 
tightest-bound comets survive these encounters over long times 
because stars need to come ever closer to destabilise the close-in 
comets. In addition, as the Sun spends more time in the cluster, 
the tidal forces from the gas and stars have more time to torque 
the comets' perihelia away from the planets. Since this time 
scale is t q cc p g a~ 2 (Duncan et al., 1987), spending more time 
in the cluster will torque comets at smaller semi-major axis. 
The same argument holds for the torquing of the comets' per- 
ihelia by stellar encounters. Secondly, ignoring the red curve 
in the upper-right panel of Fig. |4] where the Sun suffered a 
very close encounter with another star early on, the range of 
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Figure 3: A continuation of Fig. [2] The gas begins to decay after 4 Myr. 
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Figure 4: Cumulative semi-major axis for Oort clouds for various Hernquist 



clusters. Red line: t g = 2 Myr, r, 
t g = 2 Myr, r, = 6r c . Magenta: t g = 



= 3r c . Green: U = 4 Myr, r, = 3r c . Blue: 
4 Myr, r, = 6r c . Data is for sfe of 10%. 



the curves of the same colour are similar among all the pan- 
els of each model. However, when comparing the Hernquist 
and Plummer clouds, for the Hernquist model the range of the 
clouds is from approximately 200 AU to 200 000 AU (where we 
truncated it), while for Plummer the clouds are farther away, 
starting at about 800 AU. This makes sense: the central den- 
sity in the Hernquist model is higher than that in the Plummer 
model and thus as the Sun flies through the centre of the cluster 
the torquing by the gas and the encounters with the other stars 
in the Hernquist model are more violent than in the Plummer 
model. A third interesting feature is that the cumulative distri- 
butions have a similar shape and appear just different in their 
median values. It appears as if both stochastic effects from one 
simulation to the next and the number of stars in the cluster af- 
fects the final distribution. Kaib & Quinn (2008) also reported 
that their results suffered from strong stochastic effects from 
one simulation to the next. We examined the cumulative semi- 
major axis distributions for clusters with an sfe of 0.25 and con- 
clude that they are very similar to the figures shown above. 

In order to characterise the Oort clouds better, in Fig. [6] 
we have plotted a few sample Oort clouds in semi-major axis- 
pericentre space. The data are accumulated over a series of 
clusters with different N. The top two panels pertain to Hern- 
quist models while the bottom two panels are for Plummer clus- 
ters. Note that the Plummer clusters are less centrally concen- 
trated than the Hernquist clusters. Note also that unlike the tra- 
ditional Oort cloud, which has an inner edge at approximately 
2000 AU (e.g. Dones et al., 2004, Kaib & Quinn, 2008) the 
Oort cloud formed in the cluster environment can extend all the 
way down to ~200 AU so that the cluster environment produces 
an Oort cloud that is much more centrally condensed. 

The differences among the cumulative distributions are a 
measure of the concentration of the cloud. Unfortunately there 
is no definite way to define this from the cumulative distribu- 
tions of the semi-major axis. However, we might turn to the 
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Figure 5: Same as Fig.|4]but now for Plummer clusters. Once again the sfe is 
10%. 



theory of star clusters for guidance. In analogy with star clus- 
ters, a useful way to characterise the central concentration of 
the inner Oort cloud is by considering the number density of 
the comets as a function of semi-major axis. This is similar to 
counting the stars in a cluster as a function of their projected 
distance to the centre. The number density of the comets as a 
function of the semi-major axis, n{a), can be well approximated 
by a power law when far from the Sun, and the slope is usually 
-\ or -4 (Duncan et al., 1987; Dones et al., 2004; BDL6; Kaib 
& Quinn, 2008), while close to the Sun the density is usually 
flat (e.g. Brasser et al., 2010). This suggests a best fit through 
the number density profile of the form n{a) = «o(l + a/ao)~ 4 , 
where n = 3N c l(Ana*) (Dehnen, 1993; Tremaine et al., 1994) 
measures the central number density of the cloud, N c is the to- 
tal number of comets in the cloud and ao is a parameter that 
measures the central concentration. The lower its value, the 
more centrally condensed the cloud is, and probably the more 
centrally condensed or long-lived the Sun's birth cluster was. 
The number of comets as a function of semi-major axis is then 
N c (a) = N c (^f^) 3 , which is a reasonable approximation to 
the curves presented in Figs. [4] and [5] The 'half-mass' radius 
ri/2 = (2 1/3 - lr'ao ~ 3.85a . In addition N c {a ) = ±N C , so 
that most of the comets are in the 'halo' rather than in the core. 



We have performed a series of fits to the data to charac- 
terise the range and typical value of ao- While performing 
the fits, we encountered cases where a better expression was 
n(a) = «o(l + a 2 /aZ)~ 2 , where «o = N c /(n 2 a^) (Elson et al., 
1987) and N c (a) = (2N C /n) tan" 1 (a /a ) - {2N c I n){a I a ){\ + 
a 2 /ag) _1 . This distribution turns quicker from a to a~ 4 around 
«o than the Dehnen distribution. However, it is slightly more 
compressed, because the half-mass radius is at ri/2 ~ 2.26flo, 
though N c (a) = 0.182 N c . As an example, for the red curve in 
the top-right panel in Fig. |4] in which the inner Oort cloud is 
rather concentrated towards the centre, we have ao = 445 AU, 
but for the red curve in the upper-left panel ao = 825 AU. For 
the same curves using Plummer profiles, ao = 3 631 AU and 
ao = 3 825 AU respectively. The density profiles for these ex- 
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Figure 6: A few sample Oort clouds in semi-major axis-pericentre space. The 
top panels are for Hernquist clusters, the bottom two are for Plummer models. 
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Figure 7: Density profiles corresponding to the red lines in upper panels of 
Figs [4] and [5] The bullets show the data, the lines are the best fits of the form 
n(a) = «o(l + a/a y*. 



Figure 8: Top: Average value of the scaling parameter ag as a function of cluster 
membership N. Bottom: Location of Sedna in the cumulative semi-major axis 
distribution. The lines show the best linear fit through the data. The sfe is 10%. 
Error bars denote the maximum and minimum values for each quantity as a 
function of N. Left column: Hernquist. Right column: Plummer. 
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Figure 9: Same as Fig.[8]but now the sfe is 25%. 



amples, and their best fits, are given in Fig. [7] We find large 
variations for ao as a function of N for one set of values of t g 
and c. The value of ciq depends more sensitively on c than on t g , 
and decreases as either of these increase. We found no correla- 
tion between the value of ao and Q, suggesting that the previous 
two parameters have a stronger effect on the final structure of 
the inner Oort cloud. For reference, for the classical Oort cloud 
that was formed in the current Galactic environment (Dones et 
al., 2004; Kaib & Quinn, 2008) a ~ 5000 AU. 

We have plotted the values of ao, averaged over the various 
combinations of t g and c as a function of log N in the top two 
panels of Fig.|8]for an sfe of 10%, and in Fig.[9]for clusters with 
sfe of 25%. Error bars denote the maximum and minimum val- 
ues that we obtained from our data for each value of N and are 
a proxy for the values of t g , c (and Q) and how stochastic the 
dynamics is. The error bars indicate that our previous assess- 
ment of the variations the cumulative semi-major axis distribu- 
tion being due to stochastic effects was essentially correct. The 



left panels refer to Hernquist clusters while the right panels re- 
fer to Plummer models. As one can see, there is a trend for ao to 
increase with N. Best fits through the data yield a steeper slope 
for Plummer than for Hernquist in both cases. Even though the 
fitting errors are large, there is a trend for ao to increase with N. 
This trend appears to be systematic and we get back to it later. 

4.2. Sedna 

In BDL6 we concluded that the median distance of the in- 
ner Oort cloud to the Sun scales with the cluster density as 
a$o (p) 1//2 - In BDL6 we said that this combination of a and 
p is found in the torquing time of the pericentre, t q , and thus if 
the product ap 1 ^ 2 is constant, then t q must be a constant. From 
BDL6 we find for the inner edge a,- = 7 7OO(p o /lOOM pc- 3 r 1/2 
AU. Thus, to get (90377) Sedna, an inner Oort cloud dwarf 
planet with a ~ 500 AU (Brown et al., 2004), we need a density 
ofp ~ 20000M o pc- 3 . Do we see this occur in our clusters and 
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does it occur for long enough to torque objects onto Sedna-like 
orbits? 

In the bottom panels of Fig. [8] we plot the value of the cu- 
mulative semi-major axis distribution at Sedna's location (a = 
500 AU), f s , as a function of logN. The bottom-left panel 
once again refers to Hernquist clusters while the bottom-right 
panel are the Plummer models. The data is averaged over t g , 
c and Q, which are incorporated into the error bars. Once 
again these denote the maximum and minimum values that we 
observed. For the Hernquist case the data can be fit as f s = 
17.81% - 4.92% log N, with errors of 17% and 25% on the co- 
efficients. For Plummer and s — 0.1 a Sedna is found in only 
70% of the cases, but then only typically at the 1% level of the 
cloud, with no trend in N. When averaging over all N, the aver- 
age values of f s for Hernquist are 4.8% for (f g , c) = (2, 3), 7.7% 
for (2,6), 2.7% for (4,3) and 7.6% for (4,6), so that f s depends 
more sensitively on c than on t g . We report that in most cases 
f s increases as either c or t g increases, but we found no system- 
atic correlation between f s and Q. 

The trend of f s decreasing and ao increasing with N can 
be explained as follows. In order to obtain a Sedna, the Sun 
needs to (i) temporarily pass through an environment with very 
high density, (ii) pass through that part of the cluster where 
the torquing on the comet is at a maximum, or (iii) experi- 
ence a close stellar passage. From BDL6 we know that in a 
Plummer cluster the torquing on the comet is a maximum when 
r ~ 0.8ap. Unfortunately, this does not coincide with the max- 
imum density, and indeed the torque on a comet vanishes at 
the centre of the cluster. Thus, to produce Sedna in a Plum- 
mer cluster requires a close stellar passage. In contrast, for the 
Hernquist model both the density and the torque are a maxi- 
mum at the centre of the cluster, so that the only requirement 
is that the Sun passes close to the centre. However, the possi- 
ble trajectories that the Sun can have to get close to the centre 
and experience the spike in density and torque that are needed, 
decrease with increasing N because the size of the cluster itself 
increases. In other words, the orbit of the Sun needs to be more 
and more radial with increasing N and thus the initial condi- 
tions for these orbits occupy a smaller region of phase space for 
the larger clusters than for smaller ones. The figures above in- 
dicate that the location and existence of Sedna are only mildly 
dependent on the cluster parameters. 

4.3. Inclination and perihelion distribution 

Figure \W\ shows the cumulative value of the cosine of the 
orbital inclination in the ecliptic plane for various inner Oort 
clouds from Hernquist clusters. As one can see, there is little 
variation among the distributions with either N, t g or c. All 
clouds are predominantly prograde since the median value of 
cos i is larger than 0. If the inclination distribution were isotropic, 
then a cumulative distribution of cos i would be a straight line 
from to 1 as cos / goes from 1 to -1. As one can see, this 
is not the case for most distributions and thus the inner Oort 
cloud is not (yet) isotropic. That said, the median inclination is 
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Figure 10: Cumulative inclination for Oort clouds for various Hernquist clus- 
ters. Red line: t g = 2 Myr, r, = 3r c . Green: t g = 4 Myr, r, = 3r c . Blue: 
t g = 2 Myr, r, = 6r c . Magenta: t g = 4 Myr, r, = 6r c . 

lower in the inner parts of the cloud and higher in the outermost 
parts (not shown). However, even in the innermost part (less 
than 10% in cumulative semi-major axis) the median value of 
the inclination is typically between 45° and 55°. Sedna, with 
it's orbital inclination of 12°, is mostly in the bottom 5% of the 
inclination distribution. Figure QT| plots the cumulative perihe- 
lion distribution for the same clusters. Once again the different 
curves are similar in each plot. The cumulative distribution in q 
scales fairly well with In q, so that dN/dq cc q~ l and the differ- 
ential distribution is flat in q~ l . This implies that most comets 
are found with small perihelion relative to the semi-major axis. 
This is no surprise for two reasons. First the Galactic tide has 
not had time to randomise the q distribution. Second, the inner- 
most part of the Oort cloud i.e. comets with a < 2 000 AU, are 
barely affected by the Galactic tides even on Gyr time scales 
(Dones et al., 2004; BDL6), so we expect both the perihelion 
and inclination distributions to be preserved in this region. 

When examining the inclination and perihelion distributions 
for the Plummer clusters, they turn out to be very similar to 
those for the Hernquist clusters. This is to be expected because 
these distributions have not been evolved to the present age and 
thus the Galactic tide has not had time to mix the inclinations 
and eccentricities of the comets. Regarding the cumulative ec- 
centricity distribution, which is not shown, we find the cumu- 
lative distribution scales as /(e) oc e 13 , with /3 e (2, 3). A ther- 
mal distribution has the cumulative /(e) = e 2 . The Spearman 
rank correlation coefficient between the semi-major axis and 
eccentricity distributions varies from -0.07 to 0.05, suggesting 
no trend of the eccentricity to either increase or decrease with 
increasing semi-major axis. Between the semi-major axis and 
inclination distributions, the Spearman rank correlation coef- 
ficient is approximately 0.25, suggesting a weak trend for the 
inclination to increase with increasing semi-major axis. We ver- 
ified this by noting that the inner portion of the Oort cloud has 
a lower median inclination than the outer part. 
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Fisure 11: Same as Fig.fTolbut now for the perihelion distance. 

4.4. The fossilised inner Oort cloud 

In BDL6 we argued that comets with semi-major axis a < 
2 000 AU are barely affected by the Galactic tide because the 
average rotation rate of their apses is governed by planetary per- 
turbations rather than the Galactic tide. Therefore, it is likely 
that this population has preserved its original orbital structure 
since its formation during the time the Sun was in a cluster. 
We call this population the 'fossilised inner Oort cloud', and its 
most prominent member is Sedna, with 2000 CR105 - (a, q, i) = 
(220,44,22°) - and 2004 VN U2 - (a,q,i) = (350,47,26°) - 
being potential other members. The high-inclination objects 
(136199) Eris and 2004 XR I90 (Buffy) are excluded because 
their existence can be explained by other dynamical mecha- 
nisms (Gomes, 2011). Here we are interested in the orbital 
distribution of objects in this fossilised inner Oort cloud in or- 
der to directly compare it with observations. However, before 
proceeding there are a few issues that we have to consider. 

Recently Schwamb et al. (2010) attempted to constrain 
some of the properties of the Sun's birth cluster by comparing 
recent deep-field observations with Palomar of trans-Neptunian 
objects with the simulations presented in BDL6. Schwamb 
et al. (2010) concluded that the current environment of the 
outer Solar System is incompatible with the two densest clus- 
ters from BDL6 with a confidence of 95% or better; they were 
unable to reject the next-densest clusters, with central density 
po = 10000 M e pc~ 3 , for which Sedna was always at the in- 
ner edge. This should help us put some crude constraints on 
the cluster models employed in the current study that are com- 
patible with the data from Schwamb et al. (2010) and which 
are not. Re-examining the data from BDL6 we find that for the 
10 4 M pc~ 3 clusters Sedna is located at the 3% level in cu- 
mulative semi-major axis, and the innermost object has a semi- 
major axis of 218 AU. For the denser clusters the inner edge 
is closer in. However, how can we distinguish between objects 
that were placed on their current orbits by a stellar encounter 
and by planetary actions, in particular during the late epoch 
of planetary instability where Neptune most likely temporar- 
ily resided on a highly eccentric orbit (Tsiganis et al., 2005)? 



Gomes et al. (2005) and Gomes (201 1) have demonstrated that 
the combined effect of mean-motion resonances with Neptune 
and the Kozai mechanism (Kozai, 1962) can place objects from 
the Scattered Disc on orbits with a high inclination (up to 50°) 
and large pericentre distance (more than 40 AU), mimicking a 
fossilised inner Oort cloud object that was perturbed by a star. 
In addition, as Neptune migrated and its eccentricity decreased 
through dynamical friction, some objects detached from is res- 
onances or sphere of influence and now permanently reside 
on high-perihelion, high-inclination orbits (Gomes et al., 2005; 
Gomes, 2011). However, the combined effects of mean-motion 
resonances and the Kozai mechanism breaks down when the 
mean-motion resonances with Neptune begin to overlap, which 
occurs for semi-major axes larger than 200 to 250 AU (Gomes 
et al., 2005; Lykawka & Mukai, 2007). Thus any object with 
a semi-major axis longer than this value, and a pericentre dis- 
tance farther than 38 AU - beyond which chaotic diffusion stops 
(Gladman et al., 2002) - could be a fossilised inner Oort cloud 
object. Therefore, for our purpose, we eliminate objects with 
semi-major axis shorter than 250 AU and perihelia lower than 
38 AU. 

The cumulative distributions of the inclination and pericen- 
tre distance for the fossilised inner Oort cloud are very similar 
to those depicted in Figs.llOland[TThbove, apart that the distri- 
bution in q rolls over earlier than above. We decided instead to 
combine all the data from the Plummer simulations with sfe of 
10% into one plot. In Fig. [12] we plot, in the top-left panel, a his- 
togram of the semi-major axis distribution for all the objects in 
the fossilised inner Oort cloud. The stochastic behaviour of the 
semi-major axis distribution is caused by the chaotic nature of 
the simulations and the peaks and troughs should not be consid- 
ered as definitive, although the increasing trend is a systematic 
effect. The top-right panel depicts a histogram of the q distri- 
bution, and matches the q^ 1 profile discussed earlier i.e. most 
objects are trapped with a low value of q, with the median value 
being ~ 150 AU, but ~ 200 AU for the Hernquist clusters. The 
bottom-left panel displays the inclination distribution clearly 
indicating the prograde bias, while the bottom-right panel plots 
the semi-major axis vs pericentre distance. 

4.5. Efficiency 

Here we report the formation efficiency of the inner Oort 
cloud from our numerical simulations. We have plotted the for- 
mation efficiency as a function of the number of stars in the 
cluster in Fig. [13] Each filled square denotes the average effi- 
ciency for a certain combination of t g and c and the error bars 
denote the maximum and minimum. The combinations of t g 
and c are, from left to right, (2,3), (2,6), (4,3) and (4,6). The left 
panels are for Hernquist clusters with an sfe of 10% (top) and 
25% (bottom). The right panels are for Plummer clusters with 
the top having an sfe of 10% and the bottom 25%. As is clear, 
the typical efficiency does not depend strongly on the proper- 
ties of the cluster but remains at approximately 1.5%. BDL6 
reported typical efficiencies of 10% although Kaib & Quinn 
(2008) obtained efficiencies of typically 3% during their open 
cluster phase, which are in better agreement with our current 
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Figure 12: Several properties of a fossilised inner Oort cloud. This plot uses 
all the data from the Plummer model with sfe of 10%. The top-left panel plots 
a histogram of the a distribution. The top-right panel shows a histogram of 
the q distribution. The bottom-left panel shows a histogram of the inclination 
distribution while the bottom-right panel depicts a vs q. All distributions are 
normalised. 



Figure 13: Formation efficiency of the Oort cloud as a function of N. Each 
filled square lists the average efficiency and the error bars denote the maximum 
and minimum. Each filled square is for a different combination of t„ and c, 
which are, from left to right, (2,3), (2,6), (4,3) and (4,6). The left panels are 
Hernquist clusters with sfe 10% (top) and 25% (bottom). The right panels are 
for Plummer clusters with the top having an sfe of 10% and the bottom 25%. 



value. In addition, the results of Kaib & Quinn (2008) do not 
strongly depend on their cluster properties either. We cannot 
pinpoint the exact reason for the significant difference in the 
trapping efficiency between our current simulations and those 
presented in BDL6 although we can present a plausible argu- 
ment. 

Kaib & Quinn (2008) argued that their resulting inner Oort 
clouds and efficiencies were very strongly dependent on the 
closest stellar passage through their Oort clouds. If this pas- 
sage happened late then there would be little material left in the 
Oort cloud or the Scattered Disc to refill the cloud (Levison et 
al., 2004). In addition, in BDL6 the orbit of the Sun remained 
fixed in the Plummer potential and were obtained from a clus- 
ter in virial equilibrium, so that most orbits had a small radial 
excursion. Since the Sun stayed near the Plummer radius the 
fluctuations in density and flux of stellar encounters that the 
Sun experienced were limited. Here this is not the case: the 
Sun's orbit is nearly radial so that it passes close to the cluster 
centre before receding back into the halo after the initial phase 
of violent relaxation. Thus the Sun does not stay close to the 
Plummer or Hernquist radius and, especially in the latter clus- 
ters, experiences a change in background density of at least an 
order of magnitude. The passages through the centre of cluster 
occur after 1 to 3 Myr, by which time the number of comets that 
are being scattered by Jupiter and Saturn has decreased by as 
much as 80%. The close encounters with massive stars near the 
centre of the cluster and the decrease in tidal radius around the 
Sun strip many of its comets when it passes through the centre 
of the cluster. By then there are not many comets left to resup- 
ply the Oort cloud and thus the corresponding efficiency is low. 
This efficiency could possibly be increased by adding Uranus 
and Neptune and having them scatter the comets in their vicin- 
ity. However their formation mechanism and time scale is not 
well understood (Goldreich et al., 2004). Thus we have rea- 



son to believe that the efficiencies that we obtained here are in 
agreement with the expected dynamics. 

Is the existence of Sedna in agreement with a typical trap- 
ping efficiency of 1 .5%? Levison et al. (2008) and Morbidelli et 
al. (2009) estimate that there were approximately 1 000 Pluto- 
sized objects in the trans-Neptunian disc. Since the mass in this 
disc is probably comparable to that in the Jupiter-Saturn region, 
we also estimate that there were some 1 000 Pluto-sized objects 
in this region. Taking the size distribution for large Kuiper belt 
objects from Bernstein et al. (2004), with a cumulative slope of 
3.5, we estimate there were approximately 4000 Sedna-sized 
objects shortly after the formation of Jupiter and Saturn. If only 
1.5% of these ended up in the Oort cloud, this implies there are 
some 60 Sedna-sized objects in the Oort cloud, and thus we can 
expect to see one of these in the innermost 2% of the cloud. 
Thus, our low formation efficiency is compatible with the de- 
tection of one Sedna-like object. However, this estimate is un- 
certain by factors of a few because of the lack of knowledge 
of the primordial number of Sedna-like objects in the Jupiter- 
Saturn region. 

In the next section we better compare the compatibility of 
our cluster models with a single detection of Sedna. 

5. Observational Constraints on the Sedna population 

In this section we compare the orbital distribution of Sedna- 
like bodies produced in the above cluster environments to the 
observational constraints the wide-field survey of Schwamb et 
al. (2010) (hereafter S10 ) place on the Sedna population. S10 
searched 1 1 786 deg 2 down to a mean limiting R magnitude 
of ~21.3, within ±30° of the ecliptic. With the exception of 
Sedna, no new Sedna-like bodies with perihelia beyond 45 AU 
were found despite a sensitivity out to distances of ~1000 AU. 
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For further details about the survey, we refer the reader to S10. 
Here we apply the same technique and survey simulator, de- 
veloped in S10, to test whether the new cluster orbital distribu- 
tions can serve as the source of the Sedna population and are 
consistent with the single re-detection of Sedna. For each clus- 
ter environment, we compare the orbital distributions of single 
detections produced by the survey simulator to Sedna. We em- 
ploy a modified 3-dimensional Kolmogorov-Smirnov (3-D KS) 
test, detailed in S10, which simultaneously compares the semi- 
major axis, inclination, and eccentricity (a, e, i) distributions 
of single detections produced by each cluster environment to 
Sedna's orbital parameters (a=519 AU, e=0.853, /=11.9°). 

From each cluster-created orbital distribution, we randomly 
generate 3 million orbits from the final inner Oort clouds from 
our simulations, obtaining the semi-major axis, inclination, and 
eccentricity for each inner Oort cloud candidate object. The 
S10 observations probe the present-day inner Oort cloud pop- 
ulation, in particular the Sedna region, but the orbital distribu- 
tions presented here are valid only shortly after the cluster dissi- 
pates and the Sun exits the cluster. Although Sedna's orbit and 
other objects in the fossilised inner Oort cloud are dynamically 
protected from the effects of passing stars and galactic tides in 
the current solar environment, we must account for the 4.6 Gyr 
of subsequent evolution the Solar System has undergone and 
account for those effects that may have sculpted other objects 
in the Sedna region since emplacement. To avoid confusion 
with Scattered Disc and detached Kuiper belt objects, we con- 
servatively exclude objects with perihelia smaller than 50 AU 
and semi-major axes shorter than 250 AU because these could 
have dynamically interacted with Neptune. We also exclude or- 
bits with semi-major axes greater than 3 000 AU because these 
objects are capable of becoming long period comets (Kaib & 
Quinn 2009). BDL06 obtain a value of ~ 2 Gyr for the preces- 
sion frequency of Sedna; therefore we assume that the orbits of 
other Sedna-like objects have been randomised due to planetary 
effects, and thus randomly choose all other orbital angles. 

To create our sample of single detections, we randomly as- 
sign absolute magnitudes to each of our simulated inner Oort 
cloud populations. Due to the large uncertainties in the albedos 
of such a distant population, we choose to assign absolute mag- 
nitudes rather than sizes. We assume a single power-law bright- 
ness distribution where the number of objects brighter than or 
equal to a given absolute magnitude, H m . dx , is described by: 

N(H < ff ma x) = N H < Lb lO a< - H -- L6) (9) 

The brightness distribution is scaled to Nh<i.6, the number of 
bodies with an absolute magnitude brighter than or equal to 
Sedna (H = 1.6). 

For both the Hernquist and Plummer cluster models we per- 
formed the 3-D KS test for a possible range of values for a (0.2- 
0.82), including a = 0.35 and a = 0.82, the best-fit values for 
the hot (i > 5°) and cold (i < 5°) KBOs respectively (Fraser 
et al., 2010). Each single instance of the brightness distribution 



can be thought of as a separate survey, and we continued sam- 
pling the brightness distribution until we created 10000 syn- 
thetic single detections of inner Oort cloud objects for every 
cluster environment and each value of a. We found that for any 
value of a, the 3-D K-S probability of rejection is a relatively 
fiat value as Nh<\.6 increases. Therefore we restrict our analysis 
and discussion to Nh<\.6 =1. 

Table [2] lists the results from the K-S tests for each clus- 
ter environment. We have sampled the clusters according to 
sfe and number of stars, thereby crudely averaging over the 
concentration and t g . For the first column the first four letters 
refer to the potential that was used, either Hern or Plum, eOl 
or e25 refers to the sfe being either 10% of 25%, and n50 to 
nlOOO refers to the number of stars in the cluster. Except for 
the PlumeOlnlOO model, which marginally produced Sedna, 
and the Plume25nl00 model, all the orbital distributions are 
consistent with the S10 observations of a single detection of 
a Sedna-like body, and cannot be rejected at greater than the 
90% confidence level for the entire range of a that was tested. 
The majority of the simulated cluster environments reproduced 
Sedna's orbit at the inner edge of the distribution, emplacing 
the bulk of the population onto orbits with semi-major axes and 
perihelia greater than Sedna's. Unfortunately we find no direct 
correlation between cluster size and the model fit, and the ma- 
jority of the Plummer and Hernquist potential clusters are con- 
sistent with the S10 observations. Though the synthetic Hern- 
quist potential populations provide a slightly better fit to the 
data, we are unable to reject the majority of the Plummer po- 
tential clusters. The Hernquist distribution has a wider range of 
eccentricities at lower semi-major axes for the single produced 
single detections than for the Plummer model. 

Given that all models appear almost equally viable in the 
production of Sedna we must conclude that the formation of 
a Sedna population is a natural outcome of the cluster birth 
scenario independent of the initial conditions and properties 
of the cluster, and that the small variations are just the re- 
sult of stochastic effects in the cluster. Until this work the 
10 4 M pc~ 3 cluster of BDL06 was the only orbital distribu- 
tion consistent with the detection of a single Sedna-like body 
in the S10 survey. For a = 0.35 and a = 0.82 respectively, 
S10 find the best-fit values for the number of objects brighter 
than or equal to Sedna are 595^ 49 and 112^f for the 10 4 
M Q pc -3 cluster of BDL06. The Plummer and Hernquist distri- 
butions in this study resemble the 10 4 M Q pc -3 cluster orbital 
distribution, with Sedna located at the inner edge of the orbital 
distribution with many more objects having semi-major axes 
and eccentricities larger than Sedna. We expect these clusters 
would produce similar population estimates to the 10 4 M pc -3 
clusters of BDL06. To estimate the size of the Sedna popula- 
tion, we selected a few of the orbital distributions from each 
potential to examine. For each given value of a, we randomly 
assigned an absolute magnitude 50 000 times to each synthetic 
object our survey simulator created for every value of Nh<\.6- 
For each Nh<\.6 tested, accounting for survey efficiency, the 
number of simulations in which, like the real survey, one object 
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0.2 


0.35 


0.4 


0.6 


0.82 


Herne01n50 


45.85 


47.67 


47.66 


53.45 


44.91 


HerneOlnlOO 


58.28 


68.13 


66.29 


70.58 


75.1 


Herne01n250 


71.99 


72.21 


69.33 


63.66 


57.61 


Herne01n350 


67.54 


71.88 


69.73 


72.21 


64.63 


Herne01n500 


36.02 


41.74 


45.7 


57.65 


65.44 


Herne01n750 


49.45 


48.44 


47.99 


43.71 


42.34 


HerneOlnlOOO 


64.05 


62.0 


63.9 


66.9 


67.3 


Herne25n50 


42.94 


48.61 


51.57 


53.62 


52.6 


Herne25nl00 


31.35 


21.59 


17.53 


18.21 


17.69 


Herne25n250 


85.91 


86.07 


82.95 


79.91 


82.53 


Herne25n350 


55.73 


53.67 


58.71 


59.94 


75.42 


Herne25n500 


67.4 


65.2 


64.24 


61.72 


55.79 


Herne25n750 


69.63 


71.63 


68.53 


67.72 


67.73 


Herne25nl000 


77.57 


69.16 


64.24 


55.21 


38.65 


Plume01n50 


66.77 


63.63 


69.4 


69.79 


76.42 


PlumeOlnlOO 


99.26 


99.75 


99.88 


99.93 


99.99 


Plume01n250 


76.44 


72.62 


71.5 


64.69 


60.46 


Plume01n350 


82.09 


90.54 


90.09 


89.91 


91.85 


Plume01n500 


84.76 


88.47 


86.15 


82.37 


82.38 


Plume01n750 


23.69 


41.26 


41.72 


50.57 


55.14 


PlumeOlnlOOO 


82.23 


77.06 


73.69 


72.73 


61.82 


Plume25n50 


55.35 


53.9 


58.14 


58.26 


60.33 


Plume25nl00 


93.49 


96.09 


97.25 


97.87 


90.93 


Plume25n250 


63.9 


62.3 


60.99 


53.14 


48.11 


Plume25n350 


89.61 


82.36 


82.03 


77.89 


77.64 


Plume25n500 


23.45 


22.8 


20.12 


19.8 


20.2 


Plume25n750 


60.89 


47.61 


49.53 


47.61 


46.23 


Plume25nl000 


81.49 


79.0 


77.61 


67.55 


57.85 



Table 2: 3D KS test results for the cluster that produced single detections com- 
pared to Sedna's orbit. We report the confidence level we can reject the two 
distributions as drawn from the same parent population. 

on a Sedna-like orbit is detected, are tallied. Using the same 
slopes for the brightness distribution (a = 0.35 and a = 0.82), 
we find similar results than S10. In other words, a cluster envi- 
ronment could emplace on the order of a hundred to a thousand 
planetoids brighter than or equal to Sedna beyond the Kuiper 
belt, roughly consistent with the earlier estimate of approxi- 
mately 100 objects based on the efficiency of inner Oort cloud 
formation and a crude estimate as to the number of available 
objects. An order of magnitude or two more mass may reside 
in the Sedna region than exists in the present Kuiper belt. 

6. Discussion 

The fact that many clusters seem to fit the currently-known 
structure of the outer Solar System requires further elaboration. 
From the above simulations that we performed and the present 
location of Sedna, we cannot constrain the number of stars in 
the Sun's birth cluster. The only conclusion we can draw is that 
almost all the clusters are consistent with the detection of one 
Sedna-like object at the inner edge of the Oort cloud, and thus 
the production of Sedna at its current location is a generic out- 



come. The observational data marginally favours the Hernquist 
model over the Plummer clusters, even though the former pro- 
duces an inner Oort cloud that is more centrally condensed than 
the latter. In order to determine which model best corresponds 
with the reality deeper surveys, such as LSST (Large Synoptic 
Survey Telescope), which can probe the Sedna region directly, 
are needed. 

Unfortunately, our simulations remain inconclusive as to 
what value of N best fits the current observational data. Several 
previous studies favour a large value of N, typically 1 000 to 
10000, for reasons such as supernova enrichment, gas disc trun- 
cation through photo-evaporation, gas disc truncation through 
close stellar passages and the survival of the planets in a stellar 
environment (Adams, 2010). On the other hand, Gounelle & 
Meibom (2008) argue that it is unlikely that the Solar System 
formed in an environment similar to that of the Orion Nebula 
Cluster (ONC), and argued instead that the Sun is a second- 
generation star that formed after some of the stars in the first 
generation had already gone supernova. After a few Myr of 
evolution, star formation in ONC-like settings occurs mainly in 
photo-dissociation regions where the HII region created by the 
central star and the surrounding molecular gas meet, which is 
at a few parsecs from the massive star. The fate of the ONC is 
well illustrated by the 2-3 Myr old cluster NGC 2244, whose 
most massive star has the same spectral type as Q l C Ori. In 
NGC 2244, star formation is occurring in the outskirts of the 
cluster, at distances 5-10 pc from the central 06 star. Gounelle 
& Meibom (2008) argue that the Sun formed in such a sec- 
ondary environment and that some of the short-lived radio nu- 
clei that are found in meteorites were inherited from the inter- 
stellar medium rather than supernova enrichment, or by irradi- 
ation. The first generation cluster had to be large in order for 
it to have undergone several supernovae, but constraints on the 
second-generation are weaker. Unfortunately, at this stage we 
are unable to model the second-generation cluster from which 
the Sun could have formed, and instead rely on the first genera- 
tion, from which we cannot make a definitive conclusion about 
the size and membership of the Sun's cluster. We cannot rule 
out that the Sun formed in a cluster with N > 1 000 stars as 
advocated by some studies, but at the same time we cannot rule 
out clusters with N ~ 100 either. Most stars appear to form 
in clusters with N e (100, 1 000) stars (Lada & Lada, 2003; 
Adams et al., 2006), all of which are dynamically consistent 
with the observed orbital distribution of the outer Solar System. 

Given the results above the natural question to ask is whether 
or not we could have done things differently. We performed a 
comprehensive search of the literature to determine what are the 
best initial conditions and properties for the embedded clusters 
and the gas expulsion. The only thing that we can think of is 
that our models do not update the mass distribution of the gas 
as the cluster evolves, and neither does it account for a central 
cavity in the gas distribution around the heaviest star(s). Both 
of these could affect the dynamics of the cluster. To our knowl- 
edge no N-body simulations of star clusters containing gas take 
into account the motion of and the related changes in the gas. 
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Doing so would most likely require a hybrid study of N-body 
for the stars, planets and comets, and hydrodynamics for the 
gas. This would be a completely new study in star cluster dy- 
namics far beyond the scope of this paper. 

A second issue concerns the initial conditions of the cluster. 
For this study we chose to artificially have the stars reach their 
final masses at the same time, while in reality this occurs over 
a certain time span. The sequential star formation could pre- 
vent the early violent relaxation phase of the cluster and thus 
the dynamics would be milder. Adams et al. (2006) incorpo- 
rated the sequential star formation by adding the stars one by 
one to his simulations over a time span of 1 Myr. Levison et 
al. (2010) randomised the phases of the stars by integrating the 
whole cluster in the Hernquist potential for 1 Myr. But, as we 
argued earlier, there is no definitive piece of evidence that sup- 
ports either sequential star formation or having most of the stars 
in the cluster form at more or less the same time. Further study 
is needed to support either mechanism and how it would impact 
the outer parts of the solar system. 

One additional issue that we need to address is the problem 
of when is t = 0? In our simulations, t — corresponds to when 
Jupiter and Saturn have fully formed and begin to scatter the 
comets in their vicinity. Here we assume that their time scale 
of formation is short compared to the free fall time of the Sun 
in the cluster. The free fall time is given by % = {GpY 1 ^ 2 = 
1.5 (lOOM pc~ 3 /p) 1/2 Myr. A typical number density in the 
centre of the clusters is n, = 65 pc~ 3 (Carpenter, 2000), which 
implies p ~ 28 M & pc~ 3 . Taking the sfe to be 0.1 we have 
fft ~ 1 Myr, which is more or less the time when the cluster 
reaches it maximum density in Fig. [2] The formation of Jupiter 
and Saturn is thought to have taken between 1 Myr and 3 Myr 
(Lissauer & Stevenson, 2007) and observations indicate that cir- 
cumstellar gas discs have a typical lifetime of 3 Myr (Currie et 
al., 2008). Thus, it is possible that the Sun had already passed 
through the centre of the cluster at least once before the forma- 
tion of Jupiter and Saturn had been completed. Would including 
Jupiter, Saturn and the comets after this initial passage dramat- 
ically change the outcome of our simulations? We do not think 
so because the Sun will still have a radial orbit and continue to 
experience a few passages through the cluster centre before the 
gas is blown away by stellar winds. Given that the first passage 
deposits most of the comets in the cloud and is also the most 
damaging, we believe that the outcome will be similar to what 
we have presented here. 

The last issue we would like to address is whether or not our 
model predicts more Sednas. As we discussed in Section 5, the 
existence of Sedna cannot rule out the current cluster model. 
To our knowledge, the LSST will be able to place constraints 
on the model. By extrapolating expected detections from our 
models to a survey such as LSST with a limiting magnitude of 
24, and covering ±20° from the ecliptic allows us to constrain 
the slope of the size distribution of objects in the Sedna region. 
If there are more objects in this region brighter than or equal to 
Sedna, then we expect there to be of the order of several to a few 



hundred more bodies of similar size, similar to our predictions 
from Section 5. However, the final outcome depends on the 
slope of the size distribution. If the slope is steep then we expect 
there to be many more similar objects in this region, but if the 
size distribution is rather shallow then maybe only a single or 
handful of objects detectable. The next-generation surveys such 
as LSST could hopefully shed some light on this issue. 

7. Summary and conclusions 

We have performed a series of numerical simulations of the 
formation of the inner Oort cloud in an embedded cluster envi- 
ronment. The initial conditions and other relevant parameters 
of the clusters are chosen according to the most recent obser- 
vations. Specifically, we use two potential/density pairs for the 
cluster, those of Plummer (1911) and Hernquist (1990), two 
values for the central concentration of the cluster, two values 
for the star formation efficiency and the time the gas is removed, 
and a range of values of the initial subvirial kinetic energy of the 
stars. The number of stars in the clusters are 50, 100, 250, 350, 
500, 750 and 1 000. The formation of the inner Oort cloud is 
modelled by having Jupiter and Saturn scatter comets in their 
nearby vicinity. Uranus and Neptune were not included at this 
stage. We conclude the following: 

• The inner edge of the Oort cloud ranges from approxi- 
mately 100 AU to a couple hundred AU, while the outer 
edge is beyond 10000 AU. Both are virtually indepen- 
dent on the number of stars in the cluster. 

• The central concentration of the Oort cloud is measured 
by the concentration radius, ao. Half of the mass re- 
sides within ~ 3fl(). For Hernquist clusters «o ranges from 
approximately 600 AU to 1500 AU, while for Plummer 
clusters it ranges from 1500 AU to 4000 AU. This dif- 
ference is no surprise because the Hernquist clusters are 
more centrally condensed. For the current Oort cloud 
which formed in the current Galactic environment the 
concentration radius is approximately 5000 AU. 

• The concentration radius increases with increasing num- 
ber of stars in the cluster, so that more populated clusters 
form less compact inner Oort clouds. The reason for this 
is the smaller volume of phase space that the Sun can oc- 
cupy to pass close enough to the centre to experience the 
necessary high density and torquing to deposit the comets 
in the cloud. 

• The location of the dwarf planet Sedna at 500 AU is usu- 
ally within the inner 5% of the cloud for Hernquist clus- 
ters, and within the innermost 2% of the cloud for Plum- 
mer clusters, in agreement with the idea that it is at the 
inner edge of the inner Oort cloud. 

• Almost all the orbital distributions of the clusters are con- 
sistent with a single detection of a Sedna-like body, and 
cannot be rejected at greater than the 90% confidence 
level for the entire range of the size distribution slopes 



18 



that were tested. In other words, from the position of 
Sedna and a comparison with the current structure of the 
outer Solar System we cannot constrain the number of 
stars in the cluster. 

• The inner Oort clouds formed in these clusters are not 
isotropic, but show a slight prograde preference with a 
median inclination of approximately 50°. The inner part 

is significantly more prograde while the outer part is slightly 
retrograde. 

• The perihelion distribution is flat in q~ l , so that most in- 
ner Oort cloud objects have small perihelia compared to 
their semi-major axis. 

• We examine the 'fossilised inner Oort cloud', which is 
the region inside approximately 2 000 AU where the Galac- 
tic tide has barely altered the orbital elements of the comets. 
We find that the distributions in inclination, perihelion 
and semi-major axis are similar than for the whole cloud. 
In this region we expect these distributions to be pre- 
served because the Galactic tide does not randomise the 
inclinations and eccentricities. The median perihelion 
distance is ~ 150 AU for Plummer and ~ 200 AU for 
Hernquist. 

• The typical formation efficiency of the Oort cloud is 1 .5%, 
lower than previous estimates, but consistent with a sin- 
gle detection of Sedna. 
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